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ROSERO ET AL.: PARAMETER SENSITIVITY AND INTERACTIONS IN ENHANCED NO AH-LSM 

1 ABSTRACT 

2 We use sensitivity analysis to identify the parameters that are most responsible for shaping land 

3 surface model (LSM) simulations and to understand the complex interactions in three versions of 

4 the Noah LSM: the standard version (STD), a version enhanced with a simple groundwater 

5 module (GW), and version augmented by a dynamic phenology module (DV). We use warm 

6 season, high-frequency, near-surface states and turbulent fluxes collected over nine sites in the 

7 US Southern Great Plains. We quantify changes in the pattern of sensitive parameters, the 

8 amount and nature of the interaction between parameters, and the covariance structure of the 

9 distribution of behavioral parameter sets. Using Sobol’s total and first-order sensitivity indexes, 

10 we show that very few parameters directly control the variance of the model output. Significant 

11 parameter interaction occurs so that not only the optimal parameter values differ between 

12 models, but the relationships between parameters change. GW decreases parameter interaction 

13 and appears to improve model realism, especially at wetter sites. DV increases parameter 

14 interaction and decreases identifiability, implying it is overparameterized and/or 

15 underconstrained. A case study at a wet site shows GW has two functional modes: one that 

16 mimics STD and a second in which GW improves model function by decoupling direct 

17 evaporation and baseflow. Unsupervised classification of the posterior distributions of behavioral 

18 parameter sets cannot group similar sites based solely on soil or vegetation type, helping to 

19 explain why transferability between sites and models is not straightforward. This evidence 

20 suggests a priori assignment of parameters should also consider climatic differences. 

21 Key words: sensitivity analysis, land surface models, Noah LSM, model identification 

22 Index terms: 1 843 Land/atmosphere interactions, 1 805 Computational hydrology, 0550 Model 

23 verification and validation, 1873 Uncertainty assessment, 1847 Modeling 
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1 1. Introduction 

2 Like other environmental models built to support scientific reasoning and testable 

3 hypotheses to improve our understanding of the Earth system, land-surface models (LSMs) have 

4 grown in sophistication and complexity (Pitman, 2003; Niu et al., 2009). The evaluation of LSM 

5 simulations is consequently non-trivial and, especially when LSMs are to be used in predictive 

6 mode for operational forecasting, policy assessments, or decision making, demands more 

7 powerful methods for the analysis of their behavior (Saltelli, 1999; Jakeman et al., 2006; Randall 

8 et al., 2007; Gupta et al., 2008; Abramowitz et al., 2009). One such method is sensitivity analysis 

9 (SA), which is the process of investigating the role of the various assumptions, simplifications 

10 and other input (parameter) uncertainties in shaping the simulations made by a model. SA is a 

11 tool that enables exploring high-dimensional parameter spaces of complex environmental models 

12 to better understand what controls model performance (Saltelli et al., 2008). Monte Carlo-based 

13 uncertainty and SA uses the results of multiple model realizations to evaluate the range of model 

14 outcomes and identifies the input parameters that give rise to this uncertainty (Wagener et al., 

15 2001 ; Wagener and Kollat, 2007). Used to its full potential, SA weighs model adequacy and 

16 relevance, identifies critical regions in the space of the inputs, unravels parameter interactions, 

17 establishes priorities for research, and, through an interactive process of revising the model 

18 structure, leads to simplified models and increased understanding of the natural system (Saltelli 

19 et al., 2006). In this article, we inform LSM development by using sophisticated SA to guide the 

20 development of the commonly used Noah LSM (Ek et al., 2003). 

21 SA has been an underutilized resource in LSM development. Approaches to quantify 

22 ‘sensitivity’ (the rate of change in model response with respect to a factor) have very frequently 

23 been restricted to a simple exploratory analysis of the effects of factors taken one-at-a-time 
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1 (OAT), without regard for interactions between factors. Such factors can be parameters (e.g., 

2 Pitman, 1994; Gao et al., 1996; Chen and Dudhia, 2001; Trier et al., 2008), meteorological 

3 forcing, or ancillary data sets (e.g., Kato et al., 2007; Gulden et al., 2008a). OAT is illicit and 

4 unjustified for non-linear models (Saltelli, 1999; Bastidas et al., 1999; Saltelli et al., 2006). More 

5 powerful and sophisticated approaches that account, explicitly or implicitly, for parameter 

6 interactions include the factorial method (e.g., Liang and Guo, 2003; Oleson et al., 2008) and 

7 regionalized sensitivity analysis (RSA) (e.g., Bastidas et al., 2006, Demaria et al., 2007; 

8 Prihodko et al., 2008). The factorial method is a global variance-based sensitivity analysis (VSA) 

9 in which a selected number of model runs whose parameters have been perturbed by an arbitrary 

10 percent from an arbitrary reference value (default) are evaluated to identify parameters that affect 

11 the variance of the model output and to calculate a rough estimate of the low-order interactions. 

12 RSA, on the other hand, representatively samples the entire parameter space and provides a more 

13 robust assessment of the way in which parameter distributions change between subjectively 

14 defined ‘good’ and ‘bad’ (i.e., behavioral and non-behavioral) model simulations. RSA does not 

15 explicitly account for interactions between parameters and is typically applied with the sole 

16 purpose of identifying sensitive parameters to reduce the dimensionality of the calibration 

17 problem (e.g., Bastidas et al., 1 999). When RSA and VSA are used separately, both the lack of 

18 firm conclusions regarding the effect of dominant parameters (and their interactions) on the 

19 output (e.g., Bastidas et al., 2006) and the inability to draw cause-effect relationships between 

20 parameter regions and model responses (e.g., Liang and Guo, 2003) have precluded SA findings 

21 from being widely used in model development. Our approach complements the information 

22 about the posterior distributions of behavioral parameters (RSA) with information about their 

23 contributions to model variance (VSA). The Monte Carlo-based VSA that we use (the method of 
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1 Sobol') is more robust than factorial analysis because it employs a representative sample of 

2 parameter space and bypasses the perceived ‘complexities’ often associated with powerful SA 

3 methods. 

4 Our SA guides the investigation of the performance and physical realism of three 

5 versions of the Noah LSM: the standard Noah (STD), a version augmented with a simple 

6 groundwater model (Niu et al., 2007) (GW), and a version augmented with an interactive canopy 

7 model (Dickinson et al., 1998) (DV). We evaluate the performance of Noah LSM in simulating 

8 the land-surface states and fluxes at nine sites in a transition zone between wet and dry climates 

9 using the datasets of IHOP 2002 (LeMone et al., 2007). Because of the strength of the land- 

10 atmosphere coupling in transition zones (Koster et al., 2004), our work focuses on warm-season 

11 climates of the US Southern Great Plains. We document how parameter interaction and 

12 sensitivity varies with model, site, soil, vegetation, and climate. We investigate the ways in 

13 which alternate model structures affect the behavior of ‘physically meaningful’ LSM parameters, 

14 focusing on the dominant parameter interactions that dictate model response. We show that only 

15 a few parameters directly control model variance and that parameter interaction is significant. 

16 We further compare the similarity of estimated multivariate posterior distributions of behavioral 

17 parameters and their sensitivity against those obtained at other sites. We show that the changes 

18 between sites are not solely controlled by soil or vegetation types but appear to be strongly 

19 related to the climatic gradient. 

20 This paper is organized as follows. Datasets, models, and methods are described in 

21 section 2. Experimental design and driving questions are formulated in section 3. Section 4 

22 presents the patterns of sensitivity obtained via RSA and the global variance-based method of 

23 Sobol'. Section 5 presents a case study demonstrating the use of SA to understand the functional 
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1 relationships between behavioral parameters, whose interaction serves to characterize model 

2 structure and test hypotheses that regard the formulation of model. Section 6 discusses 

3 implications of the results for the transferability of parameters between locations with similar 

4 physical characteristics. Conclusions are summarized in section 7. 

5 

6 2. Models, data and methods 

7 2.1. Hydrologically enhanced versions of Noah LSM 

8 We compare the standard Noah LSM release 2.7 (STD) to a version that we equipped with a 

9 short-term phenology module (DV) and one that couples a lumped, unconfined aquifer model to 

10 the model soil column (GW). 

11 2.1.1. Noah standard release 2.7 (STD) 

12 Noah (Ek et al., 2003; Mitchell et al., 2004) is a one-dimensional, medium complexity 

13 LSM used in operational weather and climate forecasting. The model is forced by incoming 

14 short- and longwave radiation, precipitation, surface pressure, relative humidity, wind speed and 

15 air temperature. The computed state variables include soil moisture and temperature, water 

16 stored on the canopy and snow on the ground. Prognostic variables include turbulent heat fluxes, 

17 and fluxes of moisture and momentum. Noah has a single canopy layer with climatologically 

18 prescribed albedo and vegetation greenness fraction. The soil profile of Noah is partitioned into 4 

19 layers (lower boundaries at 0. 1 , 0.4, 1 .0 and 2.0 m below the surface). The vertical movement of 

20 water is governed by mass conservation and a diffusive form of the Richard’s equation. 

21 Infiltration is represented by a conceptual parameterization for the subgrid treatment of 

22 precipitation and soil moisture. Drainage at the bottom is controlled only by gravitational forces; 

23 percolation neglects hydraulic diffusivity. Direct evaporation from the top soil layer, from water 
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1 intercepted by the canopy and adjusted potential Penman-Monteith transpiration are combined to 

2 represent total evapotranspiration. The surface energy balance determines the skin temperature of 

3 the combined ground-vegetation surface. Soil-layer temperature is resolved with a Crank- 

4 Nicholson numerical scheme. Diffusion equations for the soil temperature determine ground heat 

5 flux. The Noah LSM uses soil and vegetation lookup tables for static soil and vegetation 

6 parameters such as porosity, hydraulic conductivity, minimum canopy resistance, roughness 

7 length, leaf area index, etc. 

8 2.1.2. Noah augmented with a short-term dynamic phenology module (DV) 

9 We coupled the canopy module of Dickinson et al. (1998) to STD in order to compute 

10 changes in vegetation greenness fraction that result from environmental pertubations. The 

11 module allocates carbon assimilated during photosynthesis to leaves, roots, and stems; the 

12 fraction of photosynthate allocated to each reservoir is a function of, among other things, the 

13 existing biomass density. The model also tracks growth and maintenance respiration and 

14 represents carbon storage. Unlike STD, which computes greenness fraction by linear 

15 interpolation between monthly climatological values, DV represents short-term phenological 

16 variation by allowing leaf area to vary as a function of soil moisture, soil temperature, canopy 

17 temperature, and vegetation type. DV makes vegetation fraction an exponential function of leaf 

18 area index (LAI) (Yang and Niu, 2003). Because DV links vegetation fraction to dynamic LAI, 

19 DV makes direct soil evaporation, canopy evaporation, and transpiration more responsive to 

20 environmental conditions than STD. Unlike Dickinson et al. (1998), we parameterize the effect 

21 of water stress on stomatal conductance as a function of soil moisture deficit, not as a function of 

22 soil matric potential. 

23 2.1.3. Noah augmented with a simple groundwater model (GW) 
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1 GW couples a lumped unconfined aquifer model (Niu et al., 2007) to the lower boundary 

2 of the STD soil column. In GW, water flows vertically in both directions between the aquifer and 

3 the soil column. The modeled hydraulic potential is the sum of the soil matric and gravitational 

4 potentials. The relative water head between the bottom soil layer and the water table determines 

5 either gravitational drainage or upward diffusion of water driven by capillary forces. Aquifer 

6 specific yield is used to convert the water stored in the aquifer to water table depth. When water 

7 is plentiful, the water table is within the model’s soil column; if water is insufficient to maintain 

8 a near-surface aquifer, the water table falls below the soil column. An exponential function of 

9 water table depth modifies the maximum rate of subsurface runoff (for computation of baseflow) 

10 and determines the fraction of the grid cell that is saturated at the land surface (for calculation of 

11 surface runoff) (Niu et al., 2005). The model does not explicitly consider lateral transport of 

12 groundwater between grid cells. 

13 2.2. IHOP 2002 sites and datasets 

14 We used datasets available at http://www.rap.ucar.edu/research/land/observations/ihop/ 

15 from the IHOP 2002 field campaign (Weckwerth et al. 2004; LeMone et al., 2007) to evaluate 

16 the three versions of Noah LSM at nine sites along the Kansas-Oklahoma border and in northern 

17 Texas (Fig. 1). The nine stations were sited to obtain a representative sample of the region, 

18 which spans a strong east-west rainfall gradient and Bowen ratio. We used 45 days of high- 

19 frequency, multi-sensor measurements of meteorological forcing, surface-to-atmosphere fluxes, 

20 and near-surface soil moisture and temperature. Site characteristics, soil and vegetation classes, 

21 mean meteorological values, average heat fluxes and near-surface states for the observation 

22 period are summarized in Table 1 . 

23 2.3. Model initialization and spin-up 
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1 Following Rodell et al. (2005), we initialized each of the four soil layers at 50% 

2 saturation and at the multi-annual -mean temperature. To drive the spin-up (between January 1, 

3 2000, and May 13, 2002), we used downscaled North American Land Data Assimilation System 

4 (NLDAS) (Cosgrove et al., 2003) meteorological forcing, interpolated from a 60-minute to a 30- 

5 minute time step. The models were subsequently driven by IHOP 2002 meteorological forcing 

6 from May 13, 2002, to June 25, 2002 (DOY 130 to 176). For GW, water table depth was 

7 initialized assuming equilibrium of gravitational and capillary forces in the soil profile (Niu et 

8 al., 2007). 

9 2.4. Evaluation datasets 

10 To constrain and evaluate the models, we used sensible heat flux (H), latent heat flux 

11 (LE), ground heat flux (G), ground temperature (Tg), and first layer soil moisture (SMC 5 Cm ). All 

12 data was recorded at a 30-minute time step. In situ, high-frequency flux and near-surface state 

13 measurements are an integrated response of the land surface and therefore provide useful data for 

14 examining model soundness at a specific location (Bastidas et al., 2001; Stockli et al., 2008). To 

15 score model performance, we used root mean square error (RMSE). 

16 2.5. Parameters considered in the sensitivity analysis 

17 We selected 10 soil and 10 vegetation parameters of STD that have been deemed 

18 important at similar locations (Demarty et al. 2004; Bastidas et al., 2006). We included eight 

19 parameters responsible for the phenology module and four that control the groundwater module 

20 to analyze a total of 28 and 24 parameters for DV and GW, respectively. All other coefficients in 

21 the models were kept constant at the recommended values. Default values and feasible ranges 

22 (Table 2) for all parameters were taken from the literature (e.g., Chen and Dudhia, 2001; Hogue 

23 et al., 2006). 
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2.6. Methods for sensitivity analysis 

We use the variance-based method of Sobol' (Sobol', 1993; 2001) to efficiently identify 
the factors that contribute most to the variance of a model’s response. Unlike RSA, the method 
of Sobol' deals explicitly with parameter interaction and has recently been used to quantify 
model sensitivity and parameter interactions in hydrology (e.g., Ratto et al., 2007; Tang et al., 
2006, 2007; Yatheendradas et al., 2008; van Werkhoven et al., 2008). To our knowledge, it has 
not yet been used for LSM SA. 

Since variance-based SA (VS A) does not allow a mapping of the sensitivity back to the 
parameter space, we complemented our evaluation with regionalized sensitivity analysis (RSA). 

2.6.1. Sobol' indices for global variance-based sensitivity analysis (VSA) 

Sobol' indices enable researchers to distinguish the subset of independent input factors 
X={jci, .yXj,., Xk}that account for most of the variance of the model’s response Y=f(X) either by 
themselves (first-order) or due to interaction with other parameters (higher-order). For 
completeness, here we summarize the efficient Monte Carlo-based scheme presented by Saltelli 
(2002) to compute first-order and total Sobol' sensitivity indices. 

The first-order sensitivity index (Si) represents a measure of the sensitivity of 
7 = /(x lr x 2 ,„.,x i ) (the RMSE of a model realization evaluated against observations) to 

• th 

variations in parameter x,-. 5/ is defined as the ratio of the variance of Y conditioned on the i 
factor (Vi) to the total unconditional variance (V): 

s V, V(E(Y\x,) U,-E z (y) 

‘ V(Y) V(Y) V(Y) (1) 


Where 


Uf \ .fi-^r 1 , X,. 2 5 • • • 5 ^rk )/'(u 1 5 V 2 vjXU- 1 ) , X n , xU- + i) ,...,X r £ ) 

n - 1 r= i 
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1 is obtained from products of values of f computed from the sample matrix (n model realizations 

2 long) times values of / computed from another ^-realizations matrix where all k parameters 

3 except Xi are re-sampled. 

4 The estimates of the mean squared and the total variance are computed as: 


5 


6 


E%Y)=-^f(x r ,,x r2 ,..,x rk )f(x' r] ,x' r2 ,...,x' rk ) 
n r=1 


V(y)= -i,f{x rl ,x r2 ,..,x rk ) 2 -E 2 {Y) 


( 3 ) 

( 4 ) 


7 Instead of computing all 2 k -l terms of the variance decomposition: 

s y(Y)=t I y l +'ZZ v IJ + * v ajt 

i i j>i (5) 

9 (which would require as many as n2 k model runs), in addition to estimating S), it is customary to 

10 estimate only the total sensitivity index (Sri) associated with parameter X[. Sn encompasses the 

11 effect that of all the terms in the variance decomposition that include the factor Xi have on the 

12 variance of the model’s response. S T \ is estimated by the difference between the global 

13 unconditional variance of Y and the total contribution to the variance of Y that is caused by 

14 factors other than Xi, divided by the unconditional variance: 


15 


r(y)-r(i?(y|xj tf_,.-£ 2 (y) 

v(y) v(y) 


( 6 ) 


16 


1 « J V 

U -i — 7 ^ , A X rX » X r2 »• 1 X rk \f\ X r 1 ’ X r2 ’ X ri ’ X r (i+1) ’ • • X rk ) 

n- 1 r=l 


Where 


( 7 ) 


17 is obtained from products of values of f computed from the sample matrix times the values of f 

18 computed from another matrix where only x\ is re-sampled. To spare nk simulations for each 

19 model and each site (and to show that is possible to mine regular Monte Carlo runs, customary 
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1 used in land-surface modelling), we estimated the second part of this term by interpolating the 

2 response surface Y (as explained in Section 2.7.3). 

3 A significant difference between S T \ and St points to an important role of the interactions 

4 of the i factor (at all orders) in affecting Y (Saltelli et al., 2006). Identification of such parameter 

5 interactions can help guide model development. St i are also useful to identify input factors that 

6 are non-influential, which can help reduce the dimensionality of the parameter estimation 

7 problem. If an S T \ is negligible, then it is reasonable to fix that factor to any value within its range 

8 of uncertainty, and the dimensionality of the space of input factors or model parameters can be 

9 reduced accordingly (van Werkhoven et al., 2009). 

10 2.6.2. Kolmogorov-Smirnov test for regionalized sensitivity analysis (RSA) 

11 In RSA (Homberger and Spear, 1981; Spear et al., 1994; Bastidas et al., 1999, 2006), the 

12 Kolmogorov-Smirnov (K-S) test helps identify which parameters differ when model 

13 performance is ‘behavioral’ ( B ) (‘good’, ‘acceptable’) versus ‘non-behavioral’ (B). Note that 

14 RSA does not give any information regarding the extent to which a parameter affects the 

15 variance of the model output. It instead focuses on which factors are believed to affect the 

16 occurrence of behavioral model realizations. 

17 The two-sample K-S test (two-sided version) is independently applied to the univariate 

18 marginal cumulative density functions (CDF) of each parameter to test whether the distributions 

19 of B and B are different at a particular significance level. The maximum distance between the 

20 two CDFs serves to evaluate the null hypothesis that the parameter values come from the same 

21 probability distribution. It is intuitive that, if the CDFs are dissimilar from one another (as well 

22 as different from the a priori marginal distribution of that parameter), then such parameter is 

23 ‘sensitive’ (Saltelli, 2006). 
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1 Generally, a set of subjective constraints or thresholds provides a qualitative definition 

2 for B, and model realizations are classified into either B or B. In the research presented here, we 

3 considered the calibrated parameter sets (which belong to the joint posterior probability 

4 distribution that maximizes the likelihood function) as representative of B. The B sample was 

5 defined by parameter sets whose RMSE is half a standard deviation or more away from the mode 

6 of the cost function. For robustness, for each parameter, several different samples from B were 

7 taken and compared to the one of B. We used 5% as the significance level for our RSA. 

8 2.7 Sampling strategies for sensitivity analysis 

9 We generated representative samples of model parameters using Latin Hypercube 

10 Sampling (LHS) and of the behavioral parameter sets through multi-objective calibration. We 

11 also sampled from the response surfaces obtained via LHS using inexpensive Support Vector 

12 Machines (SVM). 

13 2.7.1 Latin Hypercube Monte Carlo sampling (LHS) 

14 We ran a total of 405,000 Monte Carlo simulations sampling random parameter sets 

15 (1 5,000 samples for each model and site) to obtain a detailed representation of the range of 

16 model responses. We used LHS because it combines the strengths of stratified and random 

17 sampling to ensure that all regions of the parameter space are represented in the sample (McKay 

18 et al., 1979; Helton and Davis, 2003). LHS divides each parameter range into disjoint intervals of 

19 equal probability. From each hypercube, one sample value is randomly taken. We sampled 

20 uniformly within feasible bounds (Table 2). For each sample, we recorded the RMSE of 5 

21 criteria: H, LE, G, Tg, and SMC 5cm .To create all the matrices involved in the computation of the 

22 Si index of the Sobol' method (Eq. 1), we used a modified LHS that enables replication. We 
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1 sequentially sampled the same parameter value twice from a given hypercube, while all other 

2 parameters in the set were chosen according to the LHS. 

3 2.7.2 Multi-objective Markov Chain Monte Carlo parameter estimation technique 

4 We used the efficient Markov Chain Monte Carlo sampling strategy of Vrugt et al. 

5 (2003) to approximate the joint posterior distribution of optimal parameters. The simultaneous 

6 minimization of the RMSE of multiple criteria {H, LE, G, Tg, SMC 5cm } allowed us to constrain 

7 the models to be consistent with several types of observations and facilitated the identification of 

8 the underlying posterior distribution of physically meaningful behavioral parameter sets. It is 

9 hoped that sets from the posterior distribution cause the model to mimic the processes it was 

10 designed to represent (Gupta et al., 1 999; Bastidas et al., 2001 ; Leplastrier et al., 2002; Hogue et 

11 al., 2006). The calibration algorithm runs, in parallel, multiple chains of evolving parameter 

12 distributions to provide a robust exploration of the parameter space. These chains communicate 

13 with each other through an external population of points, which are used to continuously update 

14 the size and shape of the proposal distribution in each chain. This procedure allows an initial 

15 population of parameter sets (uniformly sampled within pre-established, feasible ranges) to 

16 converge to a stationary sample, which maximizes the likelihood function and fairly 

17 approximates the Pareto set. The Pareto set (PS) represents the multi-objective tradeoff: no 

18 member of the PS can perform better with respect to one objective without simultaneously 

19 performing worse with respect to another competing objective (Gupta et al., 1 998). We used a 

20 sample of 1 50 parameter sets to represent the posterior distribution of ‘behavioral’ parameter sets 

21 (B). 

22 2.7.3 Response surface sampling (meta-modeling) using Support Vector Machines 

23 (SVM) 
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1 To inexpensively compute the second component of the U.i term (Eq. 7) for the Sobol' 

2 method (Section 2.6.1), we approximated the response surfaces of the LSMs using SVM trained 

3 on the Monte Carlo samples. We used the SVMlight libraries of Joachims (1999) available at 

4 http://svmlight.joachims.org/. 

5 SVM are a robust learning strategy that efficiently generalizes a non-linear functional 

6 dependence between a set of multi-dimensional inputs and a vector of outputs, Y=f(X). SVM use 

7 kernels to map the complex data into high-dimensional feature spaces in which linear 

8 relationships can describe separation planes. Through a sparse set of training examples that 

9 ‘support’ the hyperplane that best fits the observed data, SVM create a model to approximate Y 

10 (Vapnik, 1 998). SVM have been successfully applied in a variety of hydrological applications 

11 (e.g., Asefa et al., 2004, 2005, 2006; Khalil et al., 2005; Kaheil et al., 2008a, 2008b). SVM 

12 perform better than artificial neural networks (Dibike et al., 2001 ; Gill et al., 2008), with lower 

13 complexity (only 3 parameters), less computational cost (faster training), and reduced sensitivity 

14 to noisy inputs and sparse training samples. 

15 Based on statistical learning theory (Smola and Scholkopf, 2004), SVM are characterized 

16 by a mechanism that avoids overfitting and that results in good generalization. A coefficient C 

17 determines the tradeoff between the complexity (or flatness) of the function / and the amount to 

18 which deviations larger than empirical error e are tolerated. In the SVM training, prediction error 

19 and model complexity are simultaneously minimized in a quadratic optimization problem, which 

20 results in a uniquely global optimum (Vapnik, 1998). We used Gaussian radial basis function 

21 kernels and determined near-optimal values of the 3 SVM hyper-parameters (the capacity C, the 

22 soft margin e and the kernel width f) using the SCE calibration algorithm (Duan et al., 1992). 


14 of 67 



ROSERO ET AL.: PARAMETER SENSITIVITY AND INTERACTIONS IN ENHANCED NO AH-LSM 

1 For detailed explanation on the estimation procedure for SVM parameters, we refer interested 

2 readers to Khalil et al. (2006). 

3 The error in the estimation of the U.i term is small. For example, an SVM (C=97.12, s= 

4 0.0025 and y= 0.0323) trained on 5,000 random samples of GW [ Y= RMSE(SMC 5 Cm )] at Site 9 

5 (and calibrated against other 5,000 testing samples) predicted the remaining independent 

6 validation samples from the Y surface with a RMSE of 0.93 [%]. The estimated variance of Y 

7 (8.857) is only 3% different than the ‘true’ variance (8.584) obtained by the Monte Carlo runs. 

8 2.8. Hierarchical clustering for comparisons of parameter distributions 

9 Unsupervised classification of behavioral parameter distributions allowed us to gain 

10 understanding about the similarities of the data across locations, specifically about the 

11 relationships between types of parameters and sites. We used clustering methods to separate the 

12 collection of marginal posterior distributions of calibrated parameters sets into groups. 

13 Agglomerative hierarchical clustering methods start with n groups (one object per group) and 

14 successively merge the two most similar groups until a single group is left. We used MATLAB’s 

15 complete linkage algorithm, in which the maximum distance between objects, one coming from 

16 each cluster, represents the smallest sphere that can enclose all objects in the two groups into a 

17 single cluster (Hair et al., 1 995). Since the distance used to measure dissimilarity between 

18 observations (e.g., Manhattan, Euclidean, etc.) may influence the membership of samples to 

19 groups, we used the cophenetic correlation coefficient to assess the quality of the linkage 

20 (Martinez and Martinez, 2002). We used dendrograms to show the links between the objects as 

21 inverted U-shaped lines, whose height represents the distance between the objects. 

22 

23 3. Driving questions and experimental design 
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1 We first ask: What are the dominant model parameters across the region? We ran a suite 

2 of Monte Carlo simulations to identify parameters that exert the greatest control on the 

3 variability of simulated fluxes and states at each IHOP site for all 3 models (STD, GW and DV). 

4 We quantify sensitivity using RSA and the method of Sobol'. Our SA guides our further 

5 investigation. 

6 We then address the question: What are the dominant interactions between model 

7 parameters, and how do these change between models? With our focus toward model 

8 development, we investigate the relationships between behavioral model parameters and quantify 

9 how they change between models using the estimates of the total-order sensitivity, the 

10 multivariate posterior parameter distributions, and the covariance structures. 

11 We finally ask: How do behavioral parameters change with dominant physical 

12 characteristics of the land? We summarize the relationships between model parameters and 

13 physical characteristics by classifying the multivariate posterior distributions according to types 

14 of soil and vegetation. Our classification provides insights into how parameters can be 

15 transferred to ungauged locations. 

16 

17 4. What parameters are sensitive? 

18 We report the results of our regionalized and variance-based sensitivity analyses, which 

19 we use in the next section to help evaluate model behavior. 

20 4.1. Regionalized sensitivity analysis (RSA) 

21 The sensitivity of a given parameter changes between models and varies by location, 

22 without an easily recognizable pattern. Table 3 presents the model parameters deemed sensitive 

23 at the 5% significance level according to a Kolmogorov-Smimov (K-S) test between samples 
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1 that drive behavioral (multi-objectively calibrated) and non-behavioral simulations. A value of 1 

2 means the parameter is sensitive, 0 means insensitive. The number of sensitive parameters is 

3 tabulated by site and by class. In all versions of Noah LSM, soil parameters tend to be more 

4 sensitive than vegetation parameters. 

5 At five of the nine sites, GW has the fewest number of sensitive parameters of any of the 

6 three models, despite having more total parameters than STD. At the dry sites (1-3), the addition 

7 of GW reduces the number of vegetation parameters that are sensitive, perhaps because when 

8 GW processes are not included, the vegetation parameters are given more weight in controlling 

9 the water balance of the soil. At site 9, despite decreasing the number of sensitive vegetation and 

10 soil parameters with respect to STD, GW does not even make the parameters associated with 

11 GW sensitive. 

12 DV most often has the highest number of sensitive parameters. The number of vegetation 

13 parameters at dry sites that are sensitive in DV is significantly higher than at wet sites. Not 

14 surprisingly, the conversion factor (gl) and the specific leaf area (sla) are sensitive everywhere. 

15 The fraction of carbon into growth (fragr ) and the water stress parameter (wtsr) are sensitive at 

16 eight sites, while the coefficient for soil respiration ( rssoil ), the optical depth ( tauhf) and the 

17 wood allocation parameter ( bf) are the least often sensitive. The variability of parameter 

18 sensitivity between even similar sites perhaps points to a lack of constraint of DV with respect to 

19 the partitioning of the carbon budget or suggests that interactions with other parameters 

20 decreases the separation of the CDFs, consequently hindering the identification of sensitive 

21 parameters (Y atheendradas et al., 2008). 

22 The widespread use of RSA in LSM development has been precluded by what has been 

23 RSA’s lack of firm conclusions because of the relatively poor agreement between the 
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1 conclusions of RSA studies. Even though we and Bastidas et al. (2006) used roughly the same 

2 version of Noah (STD) and allowed many of the same parameters to vary, we found that our 

3 RSA is at best partially consistent with the assessment of Bastidas et al. (2006). One possible 

4 reason for our divergent assessments is that we did not analyze exactly the same parameters and 

5 did not sample from exactly the same ranges as did Bastidas et al. A single parameter’s 

6 sensitivity, as defined in RSA approaches, is likely a function of the values of other model 

7 parameters; alteration of the sampled range of a single parameter can influence the sensitivity of 

8 many other parameters (Spear and Homberger, 1 980). Results presented elsewhere (Rosero et 

9 al., 2009) have shown that adjusting only a partial set of parameters alters model structure (i.e., 

10 the relationships between parameters); it follows that, with changed model structure, 

11 regionalized sensitivity would also vary. 

12 Perhaps the most important reason why information obtained through RSA has not been 

13 used to inform LSM development is that RSA provides neither information regarding the 

14 magnitude of a parameter’s influence on model output nor information about the influence that 

15 the parameter exerts through interactions. RSA only identifies a subset of parameter sets that 

16 improves (if only slightly) model output. A K-S statistic that corresponds to ‘high sensitivity’ 

17 means the parameters used to generate the behavioral model realizations are ‘a lot different’ from 

18 those used to generate the nonbehavioral realizations. ‘High sensitivity’ does not mean that a 

19 parameter is highly important to (exerts significant control over) model output. Furthermore, 

20 having invariant parameter distributions in the transition from B to B is a necessary but not 

21 sufficient condition for insensitivity (Spear and Homberger, 1980). Quantifying regionalized 

22 sensitivity has the potential to allow a mapping of the model sensitivity to parameters back to the 
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1 parameter space, but without information on the change in model’s variance, RSA appears only 

2 useful for limiting the dimensionality of the calibration problem at a specific site. 

3 RSA is distinct from the variance-based analysis of sensitivity in model output to 

4 parameter values. The method of Sobol' identifies parameters that exert significant control on the 

5 variance of model output. Sobol' indices are arguably a more useful tool for model developers 

6 and the model’s operational users. 

7 4.2 Global variance-based sensitivity analysis (VSA) 

8 Our VSA shows that there are only a few parameters that, by themselves, exert 

9 significant influence on model predictions. In contrast, parameter interaction predominates and is 

10 hence the principal mechanism for sensitivity. Figures 2, 3, and 4 present, for all sites, all 

11 considered parameters, and all models, the Sobol' first-order sensitivity indexes (Si, which is the 

12 fraction of the total variance of RMSE that can be solely attributed to the i th parameter) and the 

13 residual between Sobol' ’s total and first-order sensitivity index (Sri -St, which is the fraction of 

14 total variance that results from the interaction of the z th parameter with other parameters at all 

15 orders). The easily recognizable pattern of sensitivity can be linked to physical characteristics of 

16 the sites. 

17 4.2.1 First-order sensitivity (S,) 

18 The identified sensitive parameters (Figs. 2-4) are consistent with our expectations. When 

19 parameter influences changes as we would physically expect, we interpret these results as 

20 consistent with our hypothesis that to a first order, the models are adequately representing the 

21 site-to-site variation in these key components of the water and energy cycles. We observe that, 

22 for most sites and models, the greatest first-order control on simulated top-layer soil moisture is 

23 porosity ( maxsmc ) (Fig.4a). At dry sites 1-3, where direct evaporation is presumably a major 
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1 component of LE flux, for STD and GW, the bare soil evaporation exponent ( fxexp ) exerts the 

2 most first-order control on soil moisture. The LE flux simulated by GW at dry sites is controlled 

3 by fxexp and specific yield ( rous ), which controls depth to the water table, lai directly exerts 

4 control on transpiration and hence on the surface energy budget; at the most vegetated sites (7-9), 

5 lai consequently shapes the most the variance of H and LE for both STD and GW (Fig 2a, 3 a). 

6 The initial value lai is not important to determining DV’s simulated H and LE because, unlike 

7 the other two models, DV allows the initial value of lai to change. Instead, minimum stomatal 

8 resistance ( rcmin ) exerts the most control on DV-simulated LE flux. Two new parameters 

9 associated with DV, gl and sla, which control the calculation of lai, also exert first-order control 

10 on the simulated energy balance. In sparsely vegetated sites (1-3), the Zilintikevich coefficient 

11 (czil) plays a significant role in the variance of H. 

12 The control of parameters on model variance changes between models and between sites. 

13 The shift in dominant parameters gives clues about how the models are representing the water 

14 and energy cycles from semi-arid (1-3) to semi-humid (7-9) sites. Looking only at STD, we see 

15 that as the mean annual precipitation (MAP) at the sites increases, ficexp becomes less important 

16 to determining the top-layer soil moisture, and rejkdt, a parameter involved in determining 

17 maximum rates of infiltration, becomes more important (Fig. 4a). Expanding our view to include 

18 GW, in which surface runoff is de-emphasized and subsurface runoff is emphasized (see 

19 discussion about GW’s preferred modes of operation, Section 5), we see that although fxexp is 

20 still a key parameter exerting first-order control on SMC at the dry sites, rejkdt has little direct 

21 influence on soil moisture at the wet sites. The most sensitive parameter for SMC 5cm at sites 1-3 

22 is the aquifer’s specific yield {rous), which effectively controls whether aquifer water will be 

23 accessible to the near-surface soil. Consistent with our expectations, soil suction (psisat ), which 
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1 in GW controls upward movement of water from the aquifer to the soil, has significant control 

2 on SMC within GW but not within STD, in which psisat exerts less control over soil hydraulic 

3 behavior (Fig. 4a). 

4 As sites get progressively wetter, surface exchange coefficient czil exerts progressively 

5 less influence on simulated H, and rcmin exerts progressively more influence on H (Fig. 2a). The 

6 shift in control is consistent with our expectation that at more vegetated sites, stomatal resistance 

7 will be more important to determining the surface energy balance. As a site’s MAP increases, 

8 rcmin and lai have increasing control over simulated LE, while fxexp becomes less influential 

9 (Fig. 3a). Even at dry sites (1-3), DV tends to favor larger values of vegetation fraction (shdfac) 

10 than are prescribed by STD and GW. As a consequence, DV stands apart from GW and STD in 

11 that fxexp does not directly contribute to variance of any objective at the three driest sites (with 

12 the exception of the unvegetated, bare-soil site 1, where LE is controlled by fxexp). 

13 Examinations of S'* that are not in line with expectations may be used to help modelers 

14 diagnose problems with forcing data and/or model structure. For instance, fxexp has the highest 

15 St of simulated H and LE at site 6. Site 6 receives 800 mm of precipitation each year. We do not 

16 expect direct evaporation to be a more significant component of the LE flux at site 6 than at 

17 climatically simliar sites 4 and 5 or at the semi-arid sites 1-3. The discrepancy implies that our 

18 conceptual understanding of the physical processes at site 6 is incorrect, that the model does not 

19 adequately represent the physical processes, and/or that our forcing and/or evaluation data are 

20 faulty. 

21 In contrast to what was observed with RSA, the VS A shows a clear pattern of sensitivity 

22 across the sites. Site-to-site variation in the most sensitive parameters is not chiefly governed by 

23 soil or vegetation type but, similar to other studies (e.g, Liang and Guo, 2003; Demaria et al. 
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1 2007; van Werkhoven et al., 2008), appears to be of secondary importance when compared to the 

2 influence of climatic gradient. 

3 4.2.2 Sensitivity through interactions (Sr-Sn) 

4 Figures 2b, 3b, and 4b show that interactions between parameters are responsible for 

5 most of the variance in the models’ predicted H, LE, and SMC. If we assume that the 

6 parameterizations are correct, then the significant parameter interaction indicates model 

7 overparameterization (Saltelli et al., 2008). Overparameterization causes increased parameter 

8 interaction (observed here) and decreases the percentage of parameters that are deemed sensitive 

9 by RSA (Bastidas et al., 2006; Yatheendradas et al., 2008) (observed in the RSA described 

10 above). VSA enables the quantification of the significance of parameter interaction and the 

11 identification of the most interactive parameters. Although parameter interaction may not be an 

12 inherently negative trait (e.g., in porous media, we expect hydraulic conductivity and porosity to 

13 be functionally related), especially when there are no known functional relationships between the 

14 physical quantities that developers believe two parameters represent, parameter interaction is 

15 likely to be indication that the model is working in ways that are not consistent with the 

16 conceptual model from which the parameterizations were built. 

17 All models exhibit the most parameter interaction at the driest sites, consistent with the 

18 findings of Liang and Guo (2003) and suggesting the need to revise the formulation of all three 

19 models for semi-arid regions (Hogue et al., 2005; Rosero and Bastidas, 2007). Especially for H 

20 and SMC, GW reduces parameter interaction at the middling moisture (4-6) and semi-humid 

21 sites (7-9) (e.g., Fig.5b). GW’s reduction of parameter interaction is evidence (although by no 

22 means conclusive) that GW is more realistic than STD at sites 4-9. This is consistent with 

23 observations on the robustness of GW (Gulden et al., 2007). Conversely, GW appears to produce 
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1 significant parameter interaction at semiarid sites 1 -3, indicating a need for an improved 

2 parameterization of groundwater processes at semi-arid sites. DV parameters are much more 

3 interactive than those of STD and GW, especially at the wettest sites when simulating LE and 

4 SMC. The increased interaction between the DV-specific parameters and the rest of the 

5 conceptually unrelated STD parameters suggests DV is not functioning as its developers 

6 intended. The significant parameter interaction is consistent with the poor robustness of DV 

7 (Rosero et al., 2009). 

8 There is significant interaction between the vegetation parameters and the sensitivity of 

9 soil parameters. These results fully support the findings of Demaria et al. (2007) and Liang and 

10 Guo (2003). The significant interplay between soil and vegetation parameters implies that 

11 parameter interaction between different model parts changes both the optimal value of a given 

12 parameter and the relationships between parameters. In the next section, we use the multivariate 

13 posterior distribution of behavioral parameters to characterize the impact of adding new modules 

14 on parameter interactions. 

15 Looked at in full, we conclude that the models best represent the surface water and 

16 energy balance at the intermediate moisture and wet sites, where parameter interaction tends, 

17 within a given model, to be lowest. Because it reduces parameter interaction, GW is most likely 

18 of any of the three models to be representing the key physical processes with the most realism. 

19 

20 5. How do sensitive parameters interact and shape model behavior? 

21 We present a case study in which sensitivity analysis (SA) links model identification and 

22 model development. SA identifies parameters whose behavior merits further examination: VSA 
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1 indicates which parameters are most influential in determining model output; RSA helps us 

2 identify modes in parameter distributions that yield distinctive model behavior. 

3 5.1 A case study at Site 7 

4 At site 7, STD, DV, and GW show nearly equivalent performance when using their 

5 behavioral parameter sets (Figs. 5, 6). Such ‘equifinality’ is well documented in the hydrologic 

6 modeling literature (e.g. Beven and Freer, 2001 ; Rosero et al., 2009). In this case, distinguishing 

7 a ‘best’ model is not trivial. It requires us not only to confront the simulations with observed 

8 behavior to test for consistency but to fully understand the underlying model structures (the 

9 relationship between parameters) that make STD, GW and DV perform equally well. We show 

10 how sensitivity analysis offers the power and the ability to discriminate between model 

11 structures that conform to our physical understanding of the systems. 

12 5.2 Focus on sensitive parameters to better understand model function 

13 We follow the impact of shifted preferred values of three ‘physically meaningful’ 

14 parameters that made considerable contributions to first-order variance: porosity ( smcmax ), the 

15 muting factor for vegetation’s effect on thermal conductivity ( sbeta ), and minimum stomatal 

16 resistance ( rcmin ). The fundamental implication of our observations is that although the different 

17 optimal values of parameters are important (as found during model identification), the change in 

18 functional relationship between the parameters (the information contained in the interactions) is 

19 most relevant for purposes of model development. 

20 Figure 7 shows the marginal CDF of the posterior multivariate distribution for selected 

21 sensitive parameters in the behavioral range at Site 7. Along with the CDFs, box plots show the 

22 median and interquartile ranges of the skewed parameter distributions. The scatterplots show the 

23 variation in RMSE of LE and SMC 5cm along the range of plausible values of the parameter. The 
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trend is emphasized by fitting the points with a polynomial of minimal complexity. Even a 
cursory analysis of Figure 7 shows that the direction of ‘sensitivity’ (understood as the rate of 
change of score with parameter value) changes between models (e.g., Fig. 7a). The simulation of 
SMC 5 Cm by STD and DV degrades as porosity increases, while GW improves. We also note that, 
along the plausible range, the response can be enhanced (Fig. 7d) or become relatively 
insensitive to changes in parameter value (Fig. 7c). The identifiability of parameters (understood 
as having a clearly defined local minimum) changes between models. For example, in DV, there 
is a clear low point of the RMSE of LE along the range of values of the maximum water-holding 
capacity of the canopy ( cmcmax ), but STD and GW have less of a preference (Fig. 7c). The 
interquartile range of rcmin of STD is smaller than that of GW or DV (Fig. 7b). We can clearly 
see that different models have distinct optimal parameter values for the same physical parameter, 
implying not only that the parameters cannot be transferred between models but that the 
relationships between them are different. 


5.2.1 The role of porosity (maxsmc) 

In all three versions of Noah, higher values of maxsmc tend to decrease direct 
evaporation from the first soil layer (E DI r). E D ir is estimated as the product of Penman’s potential 
evaporation ( ET pot ), the complement of the vegetated fraction (shdfac), and the ratio of top-layer 
volumetric soil moisture (SMCi) to maxsmc: 


E dir =ET,(\- shdfac I 


/ SMQ- SMC ; 


\fx exp 


dry 


\majcsmc- SMC dry J 


( 8 ) 


SMCdry is the lowest volumetric water content that can exist in the top soil layer, and fxexp is a 
parameter ranging from 0.2 to 4. 
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In STD and DV, the error in simulating LE tends to be relatively small when maxsmc is 
low and relatively large when maxsmc is high (Fig. 7a). However, GW shows slight 
improvement in simulating LE as maxsmc increases. The tendency of STD and DV to simulate 
LE well when maxsmc is low (and direct evaporation from the soil consequently tends to be 
high) implies that STD and DV often underestimate direct evaporation at site 7. Given the same 
maxsmc value, GW more easily simulates sufficient direct evaporation, perhaps because of 
wetter soil (Rosero et al., 2009). 

In STD and DV, maxsmc controls both surface and subsurface runoff. Hydraulic 
conductivity of the soil layer ( wend) is computed using Clapp and Homberger’s method, which 
scales saturated hydraulic conductivity ( dksat ) by wetness (SMC/ maxsmc), raised to an 
exponent determined by the Clapp and Homberger parameter ( b ): 


wend = dksat 


SMC 


n26+3 


V maxsmc ) 


( 9 ) 


Lower maxsmc yields higher wend, which means water moves through the soil more quickly. 
Note that, in the case of parameterizing subsurface runoff, hydraulic conductivity is a 
representation of the speed at which the water can move laterally through the soil. In STD and 
DV, subsurface runoff ( Runoffl ) is wend times the slope. Consequently, higher values of maxsmc 
decrease Runoff2. Higher values of maxsmc also decrease surface runoff ( Runoffl ) by increasing 
the maximum rate of infiltration. Both changes increase the wetness of the surface soil. 

GW changes the way runoff is computed; maxsmc does not control surface or subsurface 
runoff in GW, eliminating two of the three ways in which maxsmc controls soil moisture. 

Runoffl is represented as an exponential function of depth to water (Niu et al., 2007): 

Runoffl = rsbmx e~ ffr>Zwr ( i m 
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1 Where rsbmx is the maximum rate of subsurface runoff, ^is the e-folding depth of saturated 

2 hydraulic conductivity, and Z WT is a variable that represents the depth to the water table. Runoff 1 

3 is a modified version of the same function used to compute subsurface runoff (Niu et al., 2005): 

4 Runofj\ = pcpdrp*{fsatmx e ~ 05 * fff * ZwT ^ ^ ^ 

5 Where pcpdrp is the effective incident water and the second term is the fraction of unfrozen grid 

6 cell that is saturated. 

7 In STD (and DV), maxsmc couples two physically unrelated (or very weakly related) 

8 processes (direct soil evaporation and lateral surface and subsurface runoff). GW decouples these 

9 processes by eliminating the dependence of parameterized lateral runoff on maxsmc. This 

10 decoupling reduces the spurious parameter interaction of maxsmc and, within GW, virtually 

11 eliminates the tradeoff between good simulation of LE and SMC5cm. GW is, in this regard, a 

12 better model for simulating fluxes at site 7. 

13 The question remains - why does GW poorly simulate SMC when maxsmc increases? 

14 maxsmc is used to compute vertical hydraulic conductivity (using the same function described 

15 above), which GW uses to control the flow of water between the aquifer and soil down a 

16 hydraulic gradient. Higher maxsmc yields lower hydraulic conductivity, which, in addition to 

17 decreasing the transfer of water between layers within the soil column, decreases the 

18 communication between the aquifer and the soil profile (that is, it decreases the flow of water 

19 between the two, increasing the potential for water to be retained near the surface). At site 7, GW 

20 best simulates SMC when high vertical hydraulic conductivity connects the aquifer and soil. 

21 Consistent with the work of others (e.g., Demaria et al., 2007), parameter values and 

22 model sensitivity to maxsmc are not consistent between sites along a climatic gradient or even 

23 within a set of sites with similar characteristics. Conclusions about model performance are 
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therefore difficult to generalize. This lack of continuity of behavior between sites is consistent 
with at least one of the following possibilities: (1) model parameterizations do not represent key 
aspects of the system and/or (2) our multi-objective calibration provided insufficient constraint 
for the estimation of behavioral parameters. We recommend use of runoff constraints before 
drawing firm conclusions regarding the physical reality of the runoff-related processes in GW. 

5.2.2 The role of SBETA 

All three models compute ground heat flux (G) using a flux-gradient relationship: 


G = DF l 


STC 1 — 7] 


( 12 ) 


0.5* ZSOIL {l) 

In which STCi is the temperature at the center of the first soil layer (0.5*ZSOIL(i)) and Ti is the 
surface temperature. DFi is the heat conductivity of the top soil layer. 

Noah assumes that, as vegetation cover increases, heat flux into the ground decreases. It 
uses parameter sbeta and the vegetated fraction ( shdfac ) to mute the thermal conductivity of the 
top layer of the soil {DFi) by: 


DF ] = DF ] * e abetmsWac , 1 3 

At site 7, the mode of the posterior probability distribution of all three models is near the 
bound of the explored parameter range (-1) (Fig. 7d). The preference for near-bound values is 
more pronounced in DV, which at site 7 tends to have shdfac values near 1.0, which puts 
downward pressure on the value of sbeta. The skewed posterior parameter distribution for all 
suggests that an even-less-negative value of sbeta would have yielded better results at site 7. 

The assumption that vegetation necessarily decreases the thermal conductivity of the top 
layer of the soil may be incorrect. If the ‘vegetation effect’ on thermal conductivity is real, the 
model underestimates the top-layer soil thermal conductivity. At site 7 (and at several other 
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1 sites), there is a clear tradeoff between H and G that is mediated by the thermal conductivity 

2 which hints at the need for revised process understanding. 

3 When comparing site 7 simulations to those of the other two ‘wet’ sites (8 and 9), we see 

4 a roughly consistent preference for near-zero values of sbeta. At the drier sites (1-6), the model’s 

5 strong preference for near-zero values of sbeta is less obvious; however, shdfac is closer to zero 

6 at these sites, which lowers the value of the muting factor described by Eq. 13. 

7 5.2.3 The role of minimum stomatal resistance ( rcmin ) 

8 rcmin controls much of the variance in H and LE, especially at wetter sites. As rcmin 

9 increases, the ratio of actual to potential evapotranspiration decreases, rcmin has a more 

10 consistent influence on the variance of H than on LE. 

11 At site 7, all three models prefer a low rcmin (Fig. 7b), which increases LE for a given 

12 potential evapotranspiration; however, GW and DV show decreased identifiability of rcmin. The 

13 mode of the rcmin distribution is higher for GW than for STD, perhaps because GW tends to 

14 have a wetter soil and more robust LE. The spread of the posterior distribution of rcmin for DV 

15 is significantly larger than that for STD, although both models share the same mode. This 

16 decrease in identifiability of parameters functionally related to lai (as is rcmin) is consistent with 

17 the added degrees of freedom allowed by DV (DV parameters gl and sla are most important in 

18 predicting lai [Fig. 2]). Because DV simulations include a wider spread of lai states, they also 

19 have a wider spread of ‘good’ rcmin values. 

20 5.3 What changes in GW to make it work better than or as well as STD at Site 7? 

21 Figure 8 shows the bivariate posterior distribution of selected behavioral parameters of 

22 STD and GW. The response surface of RMSE SMC 5 Cm changes between STD and GW (e.g., see 

23 maxsmc vs. psisat ). For GW, the shape of the posterior distributions of soil parameters that are 
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1 shared with STD is significantly different because of interaction with the GW parameters and 

2 module. Such shifts in model function affect the model covariance structure (Table 4). 

3 After multi-objective parameter estimation, at site 7, GW functions in one of two 

4 preferred modes (Fig. 8b). In the first, slightly preferred mode (ml), the parameters work 

5 together to help GW function as the developers likely intended. Strong communication between 

6 the aquifer and the soil column is supported by relatively high values of saturated hydraulic 

7 conductivity (satdk), low values of the reciprocal of the e-folding depth of hydraulic conductivity 

8 (fff), and low porosity ( maxsmc ). A relatively low surface runoff scaling factor (fsatmx ) and 

9 relatively high subsurface runoff scaling factor (rsbmx) ensure that subsurface runoff dominates 

10 surface runoff. Mimicking nature, high soil suction (psisat ) pulls water upward. A high aquifer 

11 specific yield (rous) deepens the water table (weakening the direct influence of the saturated 

12 zone on the model soil column) and transforms more water to runoff rather than to recharge. 

13 In the second mode (m2), GW adopts parameter values that make the model work as one 

14 would expect STD to function (i.e., the model functions with parameters that render GW 

15 nonfunctional) (Fig. 8b). Relatively high values of the reciprocal of the e-folding depth of 

16 hydraulic conductivity (fff) effectively seal the bottom of the soil column, limiting 

17 communication between the aquifer and the soil column; high porosity (maxsmc) decreases the 

18 vertical conductivity, further inhibiting the already poor communication between the soil and 

19 aquifer. High porosity favors decreased direct evaporation. Surface runoff is augmented by a 

20 relatively high surface runoff scaling factor fsatmx; subsurface runoff is lessened by the 

21 relatively low subsurface runoff scaling factor (rsbmx). 

22 These alternative behaviors are a possible explanation for the issue identified by Rosero 

23 et al. (2009), who showed that despite very good performance of calibrated GW, the model 


30 of 67 



ROSERO ET AL.: PARAMETER SENSITIVITY AND INTERACTIONS IN ENHANCED NO AH-LSM 

1 suffered from low robustness (high sensitivity to unmeasurable, errant parameters). The bimodal 

2 behavior of GW highlights the need to constrain GW with runoff observations, which will help 

3 identify a more consistent, better model. 

4 5.4 What changes in DV to make it work better than or as well as STD at site 7? 

5 DV and STD differ functionally in two ways: 1) DV predicts (rather than prescribes) 

6 shdfac as a function of environmental variation in moisture and radiation availability, and 2) 

7 rather than use a temporally constant, prescribed lai, DV uses shdfac to predict lai variation 

8 using a functional relationship: 

9 lai - max xlai mm , — log 1 (1 - shdfac) (14) 

V si 

10 The vegetated fraction affects all three components of LE: vegetation shades the soil, 

11 modulating direct evaporation (Edir); vegetation retains water above the soil, contributing to 

12 evaporation from the canopy (E c ); vegetation fuels transpiration (Etransp)- Because of shdfac' s role 

13 in scaling all components of evapotranspiration, a high value of conversion parameter gl yields a 

14 regime in which E c and Etransp are strongly favored over Edir. Low values of gl fix shdfac near 

15 zero and promote a regime in which Edir is the dominant component of LE. When shdfac is near 

16 zero, both E c and Etransp are minimized. At sites with sufficient vegetation, DV enables the model 

17 to correctly give more weight to Etransp- When there is little vegetation (e.g., at sites 1-3), the 

18 coupling to DV may be failing to consider special water use features associated with the semi- 

19 arid vegetation (Unland et al., 1 996). 

20 STD keeps shdfac fixed at monthly climatological values (~0.7 at site 7) and does not 

21 relate shdfac to lai. STD, unable to change the value of shdfac to shift the balance of components 
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1 of LE, favors higher lai (which decreases stomatal resistance and increases Etransp) as means for 

2 increasing total LE. 

3 When compared to STD, DV can achieve ‘good’ model performance using a wider range 

4 of values for shdfac and lai. We see this decreased identifiability of DV parameters when 

5 comparing the bivariate posterior parameter distributions of STD to those of DV at site 7 (Fig. 

6 9). The identifiability in the response surface of RMSE LE has changed (e.g. lai vs. rcmin). The 

7 decrease in identifiability of parameters that are functionally related to shdfac and/or lai can be 

8 seen across the IHOP sites (results not shown). 

9 The interplay of the parameters of the DV module also leads to changes in parameter 

10 densities of STD and DV (Fig. 9). We see additional evidence for increased interaction between 

11 parameters in DV when we note that the models’ covariance structure has been altered (Table 5). 

12 For example, rcmin and maxsmc are positively correlated in STD, but in DV they have a very 

13 slight negative correlation. 

14 Although the increased flexibility of lai and shdfac values may improve the model’s 

15 simulation of seasonal and interannual variation in surface fluxes, over time scales examined 

16 here, DV does not appear to improve the model. The constraints imposed by the turbulent and 

17 near-surface states may be insufficient for the complexity of the model and/or DV may need to 

18 be constrained with observations of carbon fluxes and plant growth. The function of the DV 

19 module may be hindered by Noah’s lack of a separate canopy layer (Rosero et al., 2009) or the 

20 absence of a more complex Ball-Berry type of stomatal conductance formulation (Niu et al., 

21 2009). 

22 

23 6. What does sensitivity analysis tell us about transferability? 
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1 LSM parameters are frequently assumed to be physical quantities, which can be 

2 measured and which have strong relationships with physical properties of the system. These 

3 assumptions imply that optimal values should not vary between models and that relationships 

4 between model parameters should be constant between models and between sites with similar 

5 physical properties. By making sets of vegetation-related parameters functions of vegetation 

6 type, LSMs contain the implicit assumption that vegetation type solely determines the ideal 

7 values of vegetation parameters. A priori assignment of parameter values based on a site’s 

8 physical characteristics (‘parameter transferability’) depends on the above assumption. Our VSA 

9 shows that interaction between parameters contributes most to a LSM’s total variance causing 

10 optimal parameter values to change significantly between models and sites. If parameters are 

11 transferable between similar sites, then the relationships between parameters should also be 

12 transferable between similar sites. The joint multivariate posterior distribution summarizes much 

13 of the information regarding the relationships between model parameters (model structure) at a 

14 particular location given observed datasets. We use the stable posterior distributions of the 

15 behavioral parameter sets to test the assumption that parameters and parameter relationships 

16 relate directly to physical characteristics. We also evaluate the extent to which climate 

17 determines the similarity of parameters between locations. 

18 6.1 Testing parameter transferability between sites using soil and vegetation types 

19 If parameters were readily ‘transferable’ between sites solely based on the vegetation 

20 class, we would expect the distributions of the vegetation parameters at two sites with the same 

21 vegetation type but different climatic regime (e.g., sites 2 and 8) to be more similar than the 

22 distributions of the same parameters at two sites with different vegetation but similar climate 

23 (e.g., sites 2 and 1). Using evidence from the marginal posterior parameter distributions of 
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1 behavioral, sensitive, ‘physically meaningful’ vegetation parameters ( rcmin , lai, and sbeta), 

2 Figure 10 shows that this expectation does not hold. Figure 10 contrasts the vegetation 

3 parameters’ posterior distributions for sites 2 and 8 (grassland) against the posterior distributions 

4 for contiguous sites 1 and 2 (dry). The distributions of rcmin, rsmax, lai and zO are more similar 

5 between sites with similar climate than they are between sites with the same vegetation (Fig. 10a, 

6 10b). hs and cmcmax show a similar lack of transferability. Only sbeta shows ‘transferability’ 

7 (i.e., there are smaller differences between the distributions from sites with the same vegetation 

8 cover) for all models (Fig. 10c). Parameter cfactr does show transferability, but only for STD. 

9 rgl could be considered ‘transferable,’ but only for DV and GW. 

10 Using the IHOP dataset, we cannot test parameter transferability using two sites with the 

11 same soil type but different climatology. We instead tested whether the distributions of soil 

12 parameters transfer along with soil or vegetation characteristics between these sites. Only fxexp 

13 appears ‘transferable’ between dry sites (1 and 2) with the same soil type (sandy clay loam) in 

14 STD and GW (Fig. 1 lb). Some parameters have more similar distributions between sites 2 and 8 

15 than between sites 1 and 2: psisat and satdk only for GW, Clapp and Homberger’s b only for 

16 STD, maxsmc for STD and GW (Fig. 1 la), and czil for STD and DV (Fig. 1 la). Perhaps this is 

17 because, despite sharing the same soil type, site 1 received more than twice the amount of 

18 precipitation observed at site 2. We also hypothesize that transferability based on soil or 

19 vegetation type is not uniformly viable because of the interactions between soil and vegetation 

20 parameters, which condition the shape of the distributions. The interactions between the soil and 

21 vegetation parameters at grassland site 2 are significantly different than the interactions at 

22 un vegetated site 1 . 
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1 We recognize that the case studies above are by no means conclusive, but they do not 

2 support the hypothesis that parameters are transferable solely based on vegetation type. The 

3 results instead suggest that land-surface model parameters are more sensitive to climatic forcing 

4 than to a specific soil or land-cover classification. Such results provide support to similar 

5 observations made about other hydrologic models (Demaria et al., 2007; van Werkhoven et al., 

6 2008) and made about the Noah LSM and using single optimal parameter sets (Hogue et al., 

7 2005; Rosero and Bastidas, 2007) . 

8 6.2 Synthesizing sensitivity to site, soil and vegetation classes by means of clustering 

9 In order to more quantitatively synthesize knowledge gained through sensitivity analysis 

10 for use at ungauged locations, we build upon the aforementioned idea of comparing the 

11 similarity of parameter distributions across sites (Rosero and Bastidas, 2007) by complementing 

12 the approach with unsupervised, agglomerative hierarchical clustering methods. 

13 For each IHOP site, we have obtained a stable, multivariate probability distribution x of 

14 behavioral parameter sets X={xi, X 2 , ...,Xi, .. .Jtk} using multi-objective parameter estimation. 

15 The marginal probability distribution for the I th parameter is Xi- To circumvent comparing 

16 between sites against each other (two at a time, as done above), we define a triangular probability 

17 distribution D\ as a reference distribution for each parameter. Z),= 1 when the value of parameter 

18 Xi is the “default” for the site. A=0 when jcj is at either edge of the feasible range. This step 

19 allows us to introduce the assumption that the parameters relate to soil and vegetation types. 

20 For each parameter, and at each site, we quantify the closeness between the cumulative 

21 distribution of the ‘optimal’ values of x\ (i.e, the marginal Xi) and the reference D x . We use the 

22 Hausdorff norm to quantify the difference Xi - A- For each model, the matrix of ‘signatures’ of 

23 the marginal distributions of k parameters at all the n evaluation sites is: 
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S can be used to identity groups of parameters that are similar between locations or to identify 
locations where groups of parameters behave alike. We use the unsupervised, agglomerative 
hierarchical clustering algorithm (described in Section 2.8) to find these groups without making 
any further assumptions about the number of groups. 

The Noah LSM works with three classes of parameters: parameters associated with soil 
type (soil parameters), parameters associated with vegetation cover type (vegetation parameters), 
and general parameters. If the previously described assumption of parameter transferability based 
on site characteristics holds (and if IHOP vegetation classifications are correct), then, given the 
set of signature vectors created using the set of vegetation parameter distributions S(x ve g) i..n), a 
clustering procedure should be able to classify similar sites in groups that resemble the IHOP 
vegetation type groupings (Table 1). Similarly, clustering of S(x so n,i..n) would result in sites 
grouped according to the IHOP soil type classification (Table 1). 

Applying a suite of distance metrics (e.g., Manhattan, Euclidean, Cosine, etc), neither soil 
nor vegetation parameters render groups of sites that partition solely based on the expected soil 
or vegetation classifications. Figure 12a shows the classification tree (dendrogram) for STD 
using the Euclidean distance, which maximizes the cophenetic correlation coefficient of the 
linkage (also shown). None of the distance metrics allowed us to classify S(x ves ,\..n ) by location 
in a way that matched the IHOP vegetation classifications. Given the subset S(jc so ii,i composed 
of the signature vectors of the 10 soil parameters at all sites, classification of the IHOP sites 
according to soil characteristics was also not feasible (Fig. 12b). Using signature vectors for 
STD, GW, and DV, some (but not all) of the distances related sites 7, 8, and 9 as having the 
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1 same soil and same vegetation type (although, because they also share the same climate type, we 

2 are unable to definitively attribute such classification to shared vegetation type). The rest of the 

3 sites do not strongly coalesce according to physical properties. For example, the pasture sites are 

4 not distinctively together; sites 5 and 6 (wheat crops) are never classified together according to 

5 vegetation (Fig 12a). Sites 1 and 2 (sandy clay loam) and sites 4 and 5 (loam) do not cluster 

6 together using soil parameters (Fig 12b). These results support the earlier findings presented 

7 here, which suggest that interaction between soil and vegetation parameters is significant, to the 

8 point that it shapes the posterior parameter distributions. These results also suggest that soil and 

9 vegetation type are not, by themselves, good physical characteristics by which to transfer 

10 parameters. 

11 Next, to account for the interdependence between soil and vegetation parameters, we 

12 classified the entire matrix S(jc SO ii Jtveg)- If parameters can be transferred based on shared 

13 vegetation and soil type, then the clustering of the entire matrix should identify groups of sites 

14 with the same vegetation and soil type (e.g, sites 7-9). Figure 12c shows a pattern (found with 

15 several distances) that is consistent across models: sites 7-9 cluster together. Sites 7, 8, and 9 

16 have also similar climates, and the classification of the sites shows strong resemblance to the 

17 climatic gradient. Given this dataset, we cannot disprove the contention that parameters can be 

18 transferred between sites that have both the same vegetation type and the same soil type. 

19 If we instead cluster S looking for groups of parameters, the expectation is that x SO ii will 

20 behave as a whole in a similar way across sites. In other words, one can produce a map of 

21 sensitivity to characterize which parameters are most similar to their default (prior distribution) 

22 and which are not. Figure 13 shows representative groupings of the behavioral, marginal 

23 posterior distributions of STD and GW parameters at all sites. Again, using a suite of distances, 


37 of 67 



DRAFT: JOURNAL OF GEOPHYSICAL RESEARCH (2009) 


1 we were unable to identify definitive clusters of soil and vegetation parameters within the set of 

2 signature vectors S, meaning that individual parameters are not sensitive in groups that relate 

3 primarily to soil or vegetation, but to a combination of both. The new GW parameters do seem to 

4 behave in a way that is similar to other soil parameters (Fig. 1 3b), which informs the estimation 

5 of these parameters for distributed applications. 

6 We conclude that the primary site-to-site control on the parameters identified as 

7 ‘vegetation’ or ‘soil’ parameters is not a site’s type of soil or vegetation used in isolation. This is 

8 consistent with the notion that given that parameters used by LSMs must represent physical 

9 processes across multiple scales, the parameters become “effective” rather than physically 

10 derived quantities (Wagener and Gupta, 2005) and that interaction between classes of parameters 

11 are very important. Our clustering analysis suggests that climate is a major control of site-to-site 

12 variation in parameter values and supports recommendations that mean climate be considered 

13 when transferring parameter values between sites (Liang and Guo, 2003; Demaria et al., 2007; 

14 vanWerkhoven et al., 2008). 

15 

16 7. Summary and Conclusions 

17 We combine two powerful sensitivity analysis methods, regionalized sensitivity analysis 

18 (RSA) and global variance-based sensitivity analysis (VSA), to inform model identification and 

19 development. We draw conclusions regarding LSM development and model assessment 

20 practices, the functioning of three versions of the widely used Noah LSM, and the a priori 

21 assignment of parameter values. Our work yields several conclusions that can be generalized to 

22 all LSM and to other environmental models and several others that are specific to the Noah 

23 LSM. 
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1 We show that VS A complements RSA for the purposes of model development. We 

2 perform VS A and RSA on a single set of Monte Carlo runs. When used together, VS A and RSA 

3 help developers understand model behavior (i.e., how and when different modules interact, how 

4 a model varies in performance between climate regimes, how a model performs differently 

5 between sites, how different models perform at the same site, etc.). SA identifies parameters 

6 whose behavior merits further examination: VS A identifies parameters responsible for model 

7 uncertainty; RSA identifies modes in parameter distributions that yield distinctive behavior. 

8 We show that the clear patterns of parameter importance identified by VSA are consistent 

9 with site-to-site variation in climate and with model-to-model changes in physical 

10 parameterization. VSA shows that parameter interactions within models exert significant control 

11 on model variance, and we show that the interactions can be traced using RSA. Shifts in 

12 parametric control on variance and covariance hint at whether a model represents the water and 

13 energy cycles in a way that is consistent with expectations. Although the optimal value of a 

14 parameter is useful information, the change in the functional relationship between parameters is 

15 more relevant for model development and hypothesis testing. We assert that to advance model 

16 development, VSA should be used in conjunction with RSA. 

17 Transfer of parameters based solely on shared vegetation type or on shared soil type is 

18 not a viable method for a priori parameter assignment. The work presented here shows that 

19 vegetation type and soil type are not the most significant contributors to site-to-site variance in 

20 optimal parameter values. Interaction between soil and vegetation parameters is significant and 

21 varies between sites; parameter interaction at least partially explains why transfer of parameters 

22 based solely on shared vegetation or soil type does not work. The primary factor controlling site- 

23 to-site variation in parameters is likely climate, although, given the dataset used here, the 
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1 combination of a site’s vegetation and soil type or some unidentified factor cannot be ruled out 

2 as the dominant controlling factor. The lack of viability of parameter transfer based solely on 

3 soil and vegetation type is a conclusion that has significant implications for the field of regional 

4 and global land surface modeling, which depends on parameter transfer based on stand-alone 

5 vegetation type and soil type as a means for a priori parameter assignment. 

6 Looking specifically at the performance of the three versions of the Noah LSM used here 

7 (STD, GW, and DV), we make several non-site-specific conclusions regarding model behavior. 

8 All three models exhibit significant parameter interaction, indicating that the models are 

9 overparameterized and/or underconstrained. All three show the least parameter interaction at the 

10 middling-moisture and wet sites and the most parameter interaction at the three driest sites. This 

11 difference suggests a need for reformulation of Noah LSM such that semi-arid regions are more 

12 realistically represented. On the whole, GW has less parameter interaction than STD (except at 

13 dry sites), indicating that it represents land-surface system with the most realism of any of the 

14 three models. GW is also least sensitive to errant parameters at the wettest sites (where 

15 groundwater is likely the most influential). DV has much more parameter interaction than STD, 

16 which provides evidence that the model is not performing as its developers intended, does not 

17 add value to STD, and/or requires additional constraint. Specific to site 7, we make the following 

18 observations: STD and DV tend to underestimate direct evaporation from the soil; GW does not 

19 (maybe because of wet soil). The assumption that vegetation decreases the thermal conductivity 

20 of the top layer of the soil is not well supported by data (this conclusion can be roughly 

21 generalized to other sites, especially the wet sites). At site 7, GW functions in one of two modes 

22 - the slightly preferred mode works in a way that mirrors what the developers likely intended; 

23 the second mode makes GW function as one might expect STD to work. Constraining runoff 
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1 may isolate the more realistic mode. GW has less spurious parameter interaction in part because 

2 it decouples direct evaporation and subsurface runoff (which are coupled via porosity in STD 

3 and DV). This decoupling appears to make the model function more realistically, with less 

4 tradeoff between the simulation of soil moisture and LE. Adding modules (DV, GW) decreases 

5 the identifiability of minimum stomatal resistance, although all three models prefer low 

6 minimum stomatal resistance (thus increasing LE for a given set of conditions). Across several 

7 sites, DV functions in one of two modes: the first emphasizes direct soil and canopy evaporation 

8 over transpiration; the second emphasizes transpiration over direct evaporation from the soil and 

9 canopy. 

10 Our approach to sensitivity analysis complements new methods for characterizing typical 

11 modes of LSM behavior (Gulden et al., 2008b; Rosero et al., 2009) within a model diagnostic 

12 framework (Gupta et al., 2008) that helps bridge the gap between model identification and 

13 development. We encourage other modeling groups to perform similar analyses with their 

14 models as a way to ensure rapid, continued improvement of our understanding and modeling of 

15 environmental processes. 

16 
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FIGURES 


Figure 1. IHOP 2002 near-surface state and flux stations. The contours show the strong east - west mean 
annual precipitation (MAP) gradient. The nine sites were located in representative land covers (see Table 
1): six on grassland of varying thickness, two on winter wheat, one on bare ground, and one on shrubland. 
The surface temperature of the dry (MAP=550 mm), sparsely vegetated sites (1-3) is mainly linked to the 
soil moisture. In contrast, the green, lush vegetation of the wet sites (7-9) (MAP=900 mm) controls the 
surface temperature. In sites 4-6 (MAP=750 mm), a mix of winter wheat and grassland, the surface 
temperature is influenced by both soil moisture and vegetation. 
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1 Figure 2. (a) First-order Sobol' sensitivity indices for the parameters of STD, GW and DV at all sites. Si 

2 stands for the individual contribution of a parameter to the variance of the RMSE of H. (b) Difference 

3 between Sobof’s total sensitivity index and Sj. S T i -Sj is the contribution to the variance through 

4 interactions with other parameters. Parameters grouped by soil and vegetation. See Table 2 for 

5 abbreviations of parameter names. Regional sensitivity patterns from semi-arid (MAP=550 mm), sparsely 

6 vegetated sites (1-3) to semi-humid (MAP=900 mm) sites (7-9) with green, lush vegetation, are easily 

7 distinguishable. 
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1 Figure 3. Same as Figure 2 but for LE. 
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1 Figure 4. Same as Figure 2 but for SMC 5cm . 
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3 
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1 Figure 5. Tradeoff LE-SMC 5cm and cumulative distribution functions (CDF) of scores of behavioral STD, 

2 GW and DV at Site 7. (a) Scatterplot in objective function space of parameter sets that maximize the 

3 likelihood function after multi-objective calibration against (H, LE, G, Tg, SMC 5cm }. CDF of root mean 

4 squared errors (RMSE) of behavioral runs evaluated against observed (b) LE, and (c) SMC 5cm . GW (dark 

5 grey), DV (light gray) perform as good as or better than STD (black). 

6 
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1 Figure 6. Uncertainty ranges of the 150 ensemble members of behavioral runs of STD (black), GW (dark 

2 gray) and DV (light gray) at Site 7. (a) Average diurnal LE. (b) Hourly SMC 5cm . Plots show the 

3 interquartile range IQR (50% of the runs) in continuous lines and the 90% confidence bounds (5% to 95% 

4 quantile) in dashed lines. Observations are shown with symbols. The spread of LE by STD is larger than 

5 that of DV, GW. However, IQR shows a very similar LE envelope. The ensemble of GW shows wetter 

6 SMC 5cm than the rest. The IQR of STD and DV are very similar, but the 90% confidence of STD has 

7 lower spread than DV. 
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1 Figure 7. Marginal cumulative distribution functions (CDF) of the posterior distribution of selected 

2 behavioral parameter sets at Site 7. (a) Porosity [maxsmc], (b) minimum stomatal resistance [rcmin], (c) 

3 maximum water holding capacity of the canopy [cmcmax], and (d) effect of the vegetation on ground heat 

4 flux [sbeta]. Along with the CDFS, the histograms and interquartile ranges are also shown. The trend in 

5 the scatterplots of RMSE of LE and SMC 5cm is shown by fitting a minimum complexity polynomial. Note 

6 that in all subpanels GW (dark grey), DV (light gray) and STD (black) are shown. 
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1 Figure 8. Multivariate posterior distribution of the behavioral parameters of STD and GW at site 7 shown 

2 for selected parameter combinations in bivariate plots. Higher density of parameter values are indicated 

3 with increasingly redder contours. The response surface of SMC 5 Cm is shown in the back; darker regions 

4 have higher errors. The bi-modal behavior of GW is signaled by ml and m2. See text for explanation. 
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1 Figure 9. Bivariate depiction of the posterior distribution of behavioral parameters of STD and DV at Site 

2 7. Higher density of parameter values are indicated with red contours. The response surface of LE is 

3 shown in the back; darker regions have higher errors. Note the significant change in the identifiability of 

4 hs and maxsmc. 
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1 Figure 10. Comparison of the marginal posterior parameter distributions of selected, sensitive vegetation 

2 parameters: (a) rcmin, (b) lai, and (c) sbeta. The total difference between parameter distributions at sites 

3 with the same vegetation cover type (Site 2 and 8) -grassland- (continuous, bright lines) is not smaller 

4 than the difference of distributions of the same parameters between contiguous sites (Site 2 and 1) 

5 (dashed, dark lines), which share the same soil type (sandy clay loam). 
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1 Figure 11. Comparison of the marginal posterior parameter distributions of selected, sensitive soil 

2 parameters: (a) maxsmc, (b) fxexp, and (c) czil. The total difference between parameter distributions at 

3 contiguous sites with the same soil type (Site 1 and 2) -sandy clay loam- (dashed, dark lines) is only 

4 sometimes smaller than the difference of distributions of the same parameters between sites with different 

5 meteorology (Site 2 and 8) (continuous, bright lines), which share the same vegetation cover type 

6 (grassland). 
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Figure 12. Clustering of sites using: (a) only the vegetation parameters of STD, (b) only the soil 
parameters of GW, and (c) both soil and vegetation parameters of GW. The similarity between marginal 
distributions of behavioral parameters at all sites is compared using different distances. The plots report 
the distance that maximizes the cophenetic correlation coefficient of the linkage. Note that neither soil nor 
vegetation parameters render groups solely based on soil or vegetation type. The clusters of all parameters 
seem to have a strong relationship with the 3 climatic zones: (1-3) semi-arid, (4-6) middling, and (7-9) 
semi-humid. 


(a) STD: only veg parameters 
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1 Figure 13. Clustering of soil (black), vegetation (gray) and GW-only parameters for the behavioral, 

2 marginal posterior distributions of (a) STD and (b) GW at all sites. The cophenetic correlation coefficient 

3 for the complete linkage for the parameters of STD and GW is 0.87 and 0.90, respectively. GW 

4 parameters seem to behave in a similar way as the soil parameters do. 

5 
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1 TABLES 

2 Table 1 . Average meteorology, near-surface states and turbulent fluxes observed during the calibration 

3 period (13 May - 25 Jun) at the nine IHOP 2002 sites. Indices of vegetation and soil classes are in 

4 parenthesis. Rainfall is cumulative over the observation period. Dry, sparsely vegetated sites (1-3) receive 

5 almost half of the amount of mean annual precipitation (MAP) than wet sites (7-9), with lush vegetation. 

6 Mean 2-m air temperature (Ta), sensible (H), latent (LE) and ground (G) heat flux, ground temperature 

7 (Tg) and soil moisture content at 5-cm (SMC 5cm ). 

8 


Site 

1 

2 

3 

4 

5 

6 

7 

8 

9 

Lat (°N) 

36.4728 

36.6221 

36.8610 

37.3579 

37.3781 

37.3545 

37.3132 

37.4070 

37.4103 

Lon (°W) 

100.6179 

100.6270 

100.5945 

98.2447 

98.1636 

97.6533 

96.9387 

96.7656 

96.5671 

Vegetation 

bare 

grassland 

sagebrush 

pasture 

wheat 

wheat 

pasture 

grassland 

pasture 

type 

ground 

(7) 

(9) 

(7) 

(12) 

(12) 

(7) 

(7) 

(7) 


(1) 









Soil type 

sandy 

sandy 

sandy 

loam (8) 

loam (8) 

clay loam 

silty clay 

silty clay 

silty clay 


clay loam 

clay loam 

loam (4) 



(6) 

loam (2) 

loam (2) 

loam (2) 


(7) 

(7) 








Rain (mm) 

154.5 

69.1 

72.4 

164.5 

173.6 

203.6 

175.4 

296.6 

250.8 

MAP (mm) 

530 

540 

560 

740 

750 

800 

900 

880 

900 

Ta (°C) 

21.4 

21.7 

22.5 

20.7 

20.7 

21.0 

20.7 

20.1 

19.9 

H (Wm' 2 ) 

70.5 

70.7 

75.7 

43.9 

51.9 

61.4 

25.9 

17.1 

27.9 

LE (Wm' 2 ) 

65.1 

76.1 

68.2 

106.2 

111.2 

97.1 

126.4 

122.8 

115.3 

G (Wm' 2 ) 

-10.4 

-6.4 

-9.3 

-2.7 

-5.1 

-7.5 

-5.6 

-12.1 

-10.5 

Tg (°C) 

24.1 

24.1 

25.8 

23.2 

21.9 

22.9 

22.3 

22.4 

22.7 

SMC 5cm 

15.4 

18.0 

7.0 

18.0 

18.1 

19.0 

33.2 

32.8 

34.0 

(%) 
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10 
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1 Table 2. Feasible ranges of Noah-LSM parameters considered in the sensitivity analysis. 


Parameter 

Description 

units 

min 

max 

Soil parameters 




maxsmc 

Maximum volumetric soil moisture 

3 -3 

m m 

0.35 

0.55 

psisat 

Saturated soil matric potential 

m m' 1 

0.1 

0.65 

satdk 

Saturated soil hydraulic conductivity 

™ 
m s 

IE-6 

IE-5 

b 

Clapp-Homberger b parameter 

- 

4 

10 

quartz 

Quartz content 

- 

0.1 

0.82 

refdk 

Used with refkdt to compute runoff parameter kdt 


0.05 

3 

fxexp 

Bare soil evaporation exponent 

- 

0.2 

4 

refkdt 

Surface runoff parameter 


0.1 

10 

czil 

Zilintikevich parameter 

- 

0.05 

8 

csoil 

Soil heat capacity 

Jm^K' 1 

1.26 

3.5 

Vegetation parameters 




rcmin 

Minimal stomatal resistance 

s rrf 1 

40 

400 

rgl 

Radiation stress parameter used in FI term of canopy resistance 


30 

100 

hs 

Coefficient of vapor pressure deficit term F2 in canopy resistance 


36 

47 

zO 

Roughness length 

m 

0.01 

0.1 

lai 

Leaf area index 

- 

0.1 

5 

cfactr 

Exponent in canopy water evaporation function 

- 

0.4 

0.95 

cmcmax 

Maximum canopy water capacity used in canopy evaporation 

m 

0.1 

2.0 

sbeta 

Used to compute canopy effect on ground heat flux 

- 

-4 

-1 

rsmax 

Maximum stomatal resistance 

s m' 1 

2,000 

10,000 

topt 

Optimum air temperature for transpiration 

K 

293 

303 

Dynamic Phenology parameters (Noah-DV only) 




fragr 

Fraction of carbon into growth respiration 

- 

0.1 

0.5 

gl 

Conversion between greenness fraction and LAI 

- 

0.1 

1.0 

rssoil 

Soil respiration coefficient 

s' 1 xlE-6 

0.005 

0.5 

tauhf 

Average inverse optical depth for 1/e decay of light 

- 

0.1 

0.4 

bf 

Parameter for present wood allocation 


0.4 

1.3 

wstrc 

Water stress parameter 


10 

400 

xlaimin 

Minimum leaf area index 

- 

0.05 

0.5 

sla 

Specific leaf area 

- 

5 

70 

Groundwater parameters (Noah-GW only) 




rous 

Specific yield 

3 -3 

mm 

0.01 

0.5 

fff 

e-folding depth of saturated hydraulic capacity 

-i 

m 

0.5 

10 

fsatmx 

Maximum saturated fraction 

% 

0 

90 

rsbmx 

Maximum rate of subsurface runoff 

ms' 1 IE-3 

0.01 

1 
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1 Table 3. Sensitive parameters according to a Kolmogorov-Smimov test between samples that drive 

2 behavioral (multi-objectively calibrated) and non-behavioral simulations. 1 stands for sensitive, 0 for 

3 insensitive. The number of sensitive parameter is tabulated by class (soil, vegetation, and new GW or 

4 DV). See Table 1 for parameter names. No clear regional pattern of sensitivity can be readily discerned. 
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1 Table 4. Spearman rank correlation coefficients between parameter sets belonging to the behavioral set 

2 for STD (up the diagonal) and GW (below the diagonal). Note the change in the covariance structure in 

3 Fig. 8. See Table 1 for abbreviations of parameter names. 
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1 Table 5. Spearman rank correlation coefficients between parameter sets belonging to the behavioral set 

2 for STD (up the diagonal) and DV (below the diagonal). Note the change in the covariance structure in 

3 Fig 9. See Table 1 for abbreviations of parameter names. 
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1 LIST OF FIGURES 

2 FIG 1 . IHOP 2002 near-surface state and flux stations. The contours show the strong east - west mean 

3 annual precipitation (MAP) gradient. The nine sites were sited in representative land covers (see Table 1): 

4 six on grassland of varying thickness, two on winter wheat, one on bare ground, and one on shrubland. 

5 The surface temperature of the dry (MAP=550 mm), sparsely vegetated sites (1-3) is mainly linked to the 

6 soil moisture. In contrast, the green, lush vegetation of the wet sites (7-9) (MAP=900 mm) controls the 

7 surface temperature. In sites 4-6 (MAP=750 mm), a mix of winter wheat and grassland, the surface 

8 temperature is influenced by both soil moisture and vegetation. 

9 FIG 2. (a) First-order Sobol' sensitivity indices for the parameters of STD, GW and DV at all sites. Si 

10 stands for the individual contribution of a parameter to the variance of the RMSE of H. (b) Difference 

11 between Sobof’s total sensitivity index and Sj. S T i -Sj is the contribution to the variance through 

12 interactions with other parameters. Parameters grouped by soil and vegetation. Regional sensitivity 

13 patterns from semi-arid (MAP=550 mm), sparsely vegetated sites (1-3) to semi-humid (MAP=900 mm) 

14 sites (7-9) with green, lush vegetation, are easily distinguishable. 

15 FIG 3. Same as Figure 2 but for LE. 

16 FIG 4. Same as Figure 2 but for SMC 5cm . 

17 FIG 5. Tradeoff LE-SMC 5cm and cumulative distribution functions (CDF) of scores of behavioral STD, 

18 GW and DV at Site 7. (a) Scatterplot in objective function space of parameter sets that maximize the 

19 likelihood function after multi-objective calibration against (H, LE, G, Tg, SMC 5cm }. CDF of root mean 

20 squared errors (RMSE) of behavioral runs evaluated against observed (b) LE, and (c) SMC 5cm . GW (dark 

21 grey), DV (light gray) perform as good as or better than STD (black). 

22 FIG 6. Uncertainty ranges of the 150 ensemble members of behavioral runs of STD (black), GW (dark 

23 gray) and DV (light gray) at Site 7. (a) Average diurnal LE. (b) Hourly SMC 5cm . Plots show the 

24 interquartile range (50% of the runs) in continuous lines and the 90% confidence bounds (5% to 95% 

25 quantile) in dashed lines. Observations are shown with symbols. The spread of LE by STD is larger than 

26 that of DV, GW. The ensemble of GW shows wetter SMC 5cm than the rest. 

27 FIG 7. Marginal cumulative distribution functions (CDF) of the posterior distribution of selected 

28 behavioral parameter sets at Site 7. (a) Porosity [maxsmc], (b) minimum stomatal resistance [rcmin], (c) 

29 maximum water holding capacity of the canopy [cmcmax], and (d) effect of the vegetation on ground heat 

30 flux [sbeta]. Along with the CDFS, the histograms and interquatile ranges are shown. The trend in the 
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1 scatterplots of RMSE of LE and SMC 5cm is shown by fitting a minimum complexity polynomial. Note 

2 that in all subpanels GW (dark grey), DV (light gray) and STD (black) are shown. 

3 FIG 8. Multivariate posterior distribution of the behavioral parameters of STD and GW at site 7 shown 

4 for selected parameter combinations in bivariate plots. Higher density of parameter values are indicated 

5 with increasingly redder contours. The response surface of SMC 5cm is shown in the back; darker regions 

6 have higher errors. The bi-modal behavior of GW is signaled by ml and m2. See text for explanation. 

7 FIG 9. Bivariate depiction of the posterior distribution of behavioral parameters of STD and DV at Site 7. 

8 Higher density of parameter values are indicated with red contours. The response surface of LE is shown 

9 in the back; darker regions have higher errors. Note the significant change in the identifiability of hs and 

10 maxsmc. 

11 FIG 10. Comparison of the marginal posterior parameter distributions of selected, sensitive vegetation 

12 parameters: (a) rcmin, (b) lai, and (c) sbeta. The total difference between parameter distributions at sites 

13 with the same vegetation cover type (Site 2 and 8) (continuous, bright lines) is not smaller than the 

14 difference of distributions of the same parameters between contiguous sites (Site 2 and 1) (dashed, dark 

15 lines), which share the same soil type. 

16 FIG 1 1 . Comparison of the marginal posterior parameter distributions of selected, sensitive soil 

17 parameters: (a) maxsmc, (b) fxexp, and (c) czil. The total difference between parameter distributions at 

18 sites with the same soil type (Site 5 and 7) (continuous, bright lines) is in general not smaller than the 

19 difference of distributions of the same parameters between contiguous sites (Site 5 and 6) (dashed, dark 

20 lines), which share the same vegetation cover type. 

21 FIG 12. Clustering of sites using: (a) only the vegetation parameters of STD, (b) only the soil parameters 

22 of GW, and (c) both soil and vegetation parameters of GW. The similarity between marginal distributions 

23 of behavioral parameters at all sites is compared using different distances. The plots report the distance 

24 that maximizes the cophenetic correlation coefficient of the linkage. Note that neither soil nor vegetation 

25 parameters render groups solely based on soil or vegetation type. The clusters of all parameters seem to 

26 have a strong relationship with the 3 climatic zones. 

27 FIG 13. Clustering of soil (black), vegetation (gray) and GW parameters for the behavioral, marginal 

28 posterior distributions of (a) STD and (b) GW at all sites. The cophenetic correlation coefficient for the 

29 complete linkage for the parameters of STD and GW is 0.87 and 0.90, respectively. GW parameters seem 

30 to behave in a similar way as the soil parameters do. 

31 
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1 LIST OF TABLES 

2 TABLE 1 . Average meteorology, near-surface states and turbulent fluxes observed during the calibration 

3 period (13 May - 25 Jun) at the nine IHOP 2002 sites. Indices of vegetation and soil classes are in 

4 parenthesis. Rainfall is cumulative over the observation period. Dry, sparsely vegetated sites (1-3) receive 

5 almost half of the amount of mean annual precipitation (MAP) than wet sites (7-9), with lush vegetation. 

6 Mean 2-m air temperature (Ta), sensible (H), latent (LE) and ground (G) heat flux, ground temperature 

7 (Tg) and soil moisture content at 5-cm (SMC 5cm ). 

8 TABLE 2. Feasible ranges of Noah-LSM parameters considered in the sensitivity analysis. 

9 TABLE 3. Sensitive parameters according to a Kolmogorov-Smimov test between samples that drive 

10 behavioral (multi-objectively calibrated) and non-behavioral simulations. 1 stands for sensitive, 0 for 

11 insensitive. The number of sensitive parameter is tabulated by class (soil, vegetation, and new GW or 

12 DV). No clear regional pattern of sensitivity can be readily discerned. 

13 TABLE 4. Spearman rank correlation coefficients between parameter sets belonging to the behavioral set 

14 for STD (up the diagonal) and GW (below the diagonal). Note the change in the covariance structure in 

15 Fig. 8. 

16 TABLE 5. Spearman rank correlation coefficients between parameter sets belonging to the behavioral set 

17 for STD (up the diagonal) and DV (below the diagonal). Note the change in the covariance structure in 

18 Fig. 9. 
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